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Quantum field theory reconciles quantum mechanics and special relativity, 
and plays a central role in many areas of physics. We develop a quantum 
algorithm to compute relativistic scattering probabilities in a massive quan- 
tum field theory with quartic self-interactions (0^ theory) in spacetime of four 
and fewer dimensions. Its run time is polynomial in the number of particles, 
their energy, and the desired precision, and applies at both weak and strong 
coupling. In the strong-coupling and high-precision regimes, our quantum 
algorithm achieves exponential speedup over the fastest known classical algo- 
rithm. 

The question whether quantum field theories can be efficiently simulated by quantum com- 
puters was first posed by Feynman three decades ago when he introduced the notion of quantum 
computers (7). Since then, efficient quantum algorithms have been developed to simulate the dy- 
namics of quantum lattice models and quantum systems with a fixed number of particles (2-7), 
but the question about quantum field theories has remained open. 

In this paper, we show that quantum computers can efficiently calculate scattering probabil- 
ities in continuum (p'^ theory to an arbitrary degree of precision. We choose 0^ theory because it 
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is among the simplest interacting quantum field theories, and thus illustrates the essential issues 
without unnecessary complications. We simulate a process in which initially well-separated 
particles with well-defined momentum scatter off each other. 

In our algorithm, we introduce several new techniques. First, we show that the field can 
be accurately represented with finitely many qubits by discretization of space via a lattice, and 
discretization of the field value at each lattice site. Analyzing spatial discretization errors is 
highly non-trivial for quantum field theories because of renormalization. We approach this 
problem using effective field theory. Secondly, we must create the initial state for the simulation. 
We do so by developing a modified version of adiabatic state preparation suitable for preparing 
non-eigenstates, such as wavepackets. Thirdly, to improve the efficiency of simulating the time 
evolution, we show that Suzuki-Trotter formulae converge faster in cases where the underlying 
Hamiltonians have spatial locality. These techniques may be of independent interest, beyond 
their application to simulating quantum field theory. 

No previous paper has addressed the quantum computation of scattering amplitudes or the 
convergence of quantum simulations to the continuum limit of a quantum field theory. The is- 
sue of gauge symmetries in qubit representations of lattice field theories has been studied (8), 
and there is an extensive literature on how experimentally to construct Hamiltonians that ap- 
proximate lattice gauge theories, in systems of atoms or superconducting qubits (see, for exam- 
ple, (9-19)). These previous studies are on the experimental analog implementation of lattice 
Hamiltonians, whereas the present work addresses digital simulation, with explicit consider- 
ation of convergence to the continuum, and efficient preparation of wavepacket states for the 
computation of dynamical quantities such as scattering probabilities. 

The input to our algorithm is a list of the momenta of the incoming particles, and the output 
is a list of the momenta of the outgoing particles produced by the physical scattering process. At 
relativistic energies, the number of outgoing particles may differ from the number of incoming 
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particles. However, because we consider only the case of non-zero particle mass, the number 
of outgoing particles is at most linear in the center-of-mass energy of the incoming particles. 
In accordance with quantum mechanics, the incoming momenta do not uniquely determine the 
outgoing momenta, but rather a probability distribution over possible outcomes. Upon repeated 
runs, our quantum algorithm samples from this distribution. We quantify the precision of our 
simulation by demanding that the probability of a given outcome from the simulation differ 
from the true physical probability by no more than ±e. 

The scattering processes simulated closely match experiments in particle accelerators, which 
are the standard tools to probe uniquely quantum field-theoretical effects. The problem of cal- 
culating the scattering amplitudes, encoded in an object called the ^'-matrix, has consequently 
been well studied. 

In complexity theory, the efficiency of an algorithm is judged by how its computational 
demands scale with the problem size or some other quantity associated with the problem's 
intrinsic difficulty. An algorithm with polynomial-time asymptotic scaling is considered to 
be feasible, whereas one with super-polynomial (typically, exponential) scaling is considered 
infeasible. This classification has proved to be a very useful guide in practice. Our results 
can be roughly summarized as follows: the calculation of quantum field-theoretical scattering 
amplitudes at high precision or strong coupling is infeasible on classical computers with known 
techniques but feasible on quantum computers. 

Traditional calculations of QFT scattering amplitudes rely upon perturbation theory, namely, 
a series expansion in powers of the coupling (the coefficient of the interaction term), which is 
taken to be small. A powerful and intuitive way of organizing this perturbative expansion is 
through Feynman diagrams, in which the number of loops is associated with the power of the 
coupling. A reasonable measure of the computational complexity of perturbative calculations is 
therefore the number of Feynman diagrams, which is determined by combinatorics, and grows 
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factorially with the number of loops and the number of external particles. 

If the coupling constant is insufficiently small, the perturbation series does not yield correct 
results. In (f)^ theory, for D = 2,3 spacetime dimensions, by increasing the coupling Aq, one 
eventually reaches a quantum phase transition at some critical coupling Ac (20-22). In the pa- 
rameter space near this phase transition, perturbative methods become unreliable; this region is 
referred to as the strong-coupling regime. There are then no known feasible classical methods 
for calculating scattering amplitudes, although lattice field theory can be used to obtain static 
quantities, such as mass ratios. Even at weak coupling, the perturbation series is not conver- 
gent, although it is asymptotic (23-25). Including higher-order contributions beyond a certain 
point makes the approximation worse. There is thus a maximum possible precision achievable 
perturbatively. 

We find that the number of quantum gates, Gweak> needed to sample from scattering proba- 
bilities in weakly coupled, (d + 1) -dimensional 0^ theory with accuracy ±e scales as follows^: 



The asymptotic scaling of the number of gates used to simulate the strongly coupled theory 
is summarized in Table 1. 

Although quantum field theory is typically expressed in terms of Lagrangians, and within the 
interaction picture, our algorithm is more naturally described in the formalism of Hamiltonians, 
and within the Schrodinger picture. We start by defining a lattice (f)* theory, and subsequently 
address convergence to the continuum theory. (In D — A, the continuum limit is believed to 
be the free theory. Nevertheless, since the coupling shrinks only logarithmically, scattering 
processes for particles with small momenta in lattice units are still interesting to compute.) 

'/(n) = o{g{n)) if and only if lim„^oo f{n)/g{n) = 0. In the case of e scaling it is of course 1/e thatis taken 
to infinity. We have used little-o notation to simplify our exposition. For more technical detail, see Appendix A. 
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Table 1 : The asymptotic scaling of the number of quantum gates needed to simulate scattering in 
the strong-coupling regime in one and two spatial dimensions is polynomial in the momentum 
of the incoming pair of particles, Ac — Aq, the distance from the phase transition, and nout, the 
maximum kinematically allowed number of outgoing particles. The notation /(n) = 0{g{n)) 
means f{n) — 0{g{n) log'^(n)) for some constant c. 

Let Q = alA , that is, an L x . . . x L lattice in d spatial dimensions with periodic boundary 
conditions and lattice spacing a. The number of lattice sites is V = L'^. For each x e Jl, 
let 0(x) be a continuous, real degree of freedom — interpreted as the field at x — and 7r(x) 
the corresponding canonically conjugate variable. In canonical quantization, these degrees of 
freedom are promoted to Hermitian operators with the commutation relation 



As is standard in quantum field theory, we use units with h — c — 1. (f)^ theory on the lattice Q 
is defined by the Hamiltonian 



where Va0 denotes a discretized derivative, that is, a finite-difference operator. 

We represent the state of the lattice field theory by devoting one register of qubits to store 
the value of the field at each lattice point. Each 0(x) is in principle an unbounded contin- 
uous variable. To represent the field at a given site with finitely many qubits, we cut off 
the field at a maximum magnitude ^max and discretize it in increments of 5,^. This requires 
Til, — O(log(0inax/ (^(/.)) quWts per site. Note that this field discretization is a separate issue from 
the spatial discretization via the lattice VL. 
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Let IV') be any state such that (V'li? IV') < E. The probabiUty distribution over 0(x) defined 
by I'il^) (for any x G ^2) has a very low probability^ for |0(x)| to be much larger than 0{VE). 
Thus, a cutoff ^max — O {^■\J gd^i^ suffices to ensure fidelity 1 — e to the original state 1-0). 
One can prove this by bounding (^/)|0(x)|?/)) and (V'|0^(x)|z/)) as functions of E and applying 
Chebyshev's inequality (§A.4). To choose 5^, note that the eigenbasis of a'^7r(x) is the Fourier 
transform of the eigenbasis of 0(x). Hence, discretizing 0(x) in units of (5<^ is equivalent to 
introducing the cutoff — TTmax < 7r(x) < TTmax, where TTmax = -itr- By bounding the expec- 

<X 0(j) 

tations of 7r(x) and 7r^(x), one finds that it suffices to choose TTmax = O ^y^^^)' ^^'^ ^^^^ 
». = 0(log(|5)). 

We now turn to the main three tasks of quantum simulation: preparing an initial state, simu- 
lating the time evolution e~*^*, and measuring final observables. We discuss simulation of time 
evolution first, as it is used in all three tasks. The unitary operator e~*^* can be approximated 
by a quantum circuit of 0((iV)^^^/^'^) gates implementing a k^^-order Suzuki-Trotter formula 
of the type described in (4, 26). This near-linear scaling with t has long been known. The scal- 
ing with V is a consequence of the locality^ of H (§A.6) and appears not to have been noted 
previously in the quantum algorithms literature. 

To simulate scattering, one needs to prepare an initial state of particles in well-separated 
wavepackets. We do so by preparing the vacuum of the Aq = theory, exciting wavepack- 
ets, and then adiabatically turning on the coupling Aq. Let be the Hamiltonian obtained 
by setting Aq = in H. defines an exactly solvable model in which the particles are non- 
interacting. The vacuum (ground) state | vac(O)) of is a multivariate Gaussian wavefunction 
in the variables {^(x) |x e Q}, and can therefore be prepared using the method of Kitaev and 
Webb (27). The asymptotic scaling of the Kitaev-Webb method is dictated by the computa- 

^For Ao > one has a tighter bound. In this case it is unhkely for |0(x)| to be much larger than 0{E^/^) 
(§A.4). 

couples only nearest-neighbor sites, via the (Va<?i>)^ term. 
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tion of the LDL^ decomposition of the covariance matrix, which can be done classically in 

Q(y2.376) jjjj^g ^jti^ (28,29). 

In analogy with the familiar case of the harmonic oscillator, one can define creation and 
annihilation operators Qp and aj, such that if = J2per L~'^ojpalap + E^^^ 1 , where T = gZ^ 
is the momentum-space lattice corresponding to Q, Up = y^'mp + ^ Xlj=i ^i'^^ (^) ' ^^^^ 
is an irrelevant zero-point energy. The operator a)^ can be interpreted as creating a (completely 
delocalized) particle of the non-interacting theory with momentum p and energy Up. 

The (unnormalized) state 0(x)|vac(O)) is interpreted as a single particle localized at x (see, 
e.g., {30)). Because ap|vac(0)) = 0, 0(x)|vac(O)) = a](.|vac(0)), where 

The operator 

4 = ^(^)E«'^W«x (5) 

creates a wavepacket with position-space wavefunction ^. (77 (V') is a normalization constant, 
chosen so that [a^,a^] = 1.) is not unitary, so it cannot be directly implemented by a 
quantum circuit. Instead, we introduce an ancillary qubit and let 

H^ = al®\l){Q\+a^®\Q){l\. (6) 

One can verify that e~*^''''^/^|vac(0))|0) = — ia|^|vac(0))|l). Using a high-order Suzuki-Trotter 
formula (4, 26), we can construct an efficient quantum circuit approximating the unitary trans- 
formation e"*^'/-'^/^. Applied to |vac(0)), this circuit yields the desired state up to an irrelevant 
global phase and an unentangled ancillary qubit, which can be discarded. We repeat this process 
for each incoming particle desired. 

Because we wish to create localized wavepackets, we can choose V'(x) to have bounded 
support. Expanding in terms of the operators and tt yields an expression of the form 
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alp — X^xeQ [/(^)^(^) + where /(x) and g{x.) are exponentially decaying with 

characteristic length scale 1/mo outside the support of ip. Thus, and can be exponentially 
well approximated by linear combinations of the operators (f) and tt on a local region of space, 
and the complexity of simulating e"'^*'^/^ does not scale with the volume V. Furthermore, 
provided the initial wavepackets are separated by a distance that is large compared with 1/mo, 
the preparation of each additional wavepacket leaves the existing wavepackets almost perfectly 
undisturbed. 

At this point, we have finished constructing wavepackets of the non-interacting theory. We 
next use a Suzuki-Trotter formula to construct a quantum circuit simulating the unitary transfor- 
mation induced by a time-dependent Hamiltonian in which the coupling constant is gradually 
increased from zero to its final value, Aq. By the adiabatic theorem, sufficiently slow turn-on 
ensures that no stray particles are created during this process, provided particle creation costs 
energy, that is, the particles have non-zero mass. In the free theory, the particle mass is mo. In 
the interacting theory, with fixed mo and sufficiently large Aq, the mass vanishes. This marks 
the location of the ^ —0 symmetry-breaking transition. In this paper we restrict our attention 
to simulations within the symmetric phase, although we do consider systems arbitrarily close to 
the phase transition, as these should be particularly hard to simulate classically. 

As Eq. 5 shows, wavepackets are not eigenstates of H^^K During the adiabatic turn-on, the 
different eigenstates acquire different dynamical phases. Thus, as the wavepacket time evolves, 
it propagates and broadens. This behavior is undesirable in our simulation, because we do 
not wish the particles to collide and scatter before the coupling reaches its final value. We 
therefore introduce backward time evolutions governed by time-independent Hamiltonians into 
the adiabatic state-preparation process to undo the dynamical phases. Specifically, let H{s) 
parameterize the adiabatic time evolution, with H{0) = H^°^ and H{1) = H. We divide the 
adiabatic preparation into J steps, with Uj denoting the unitary time evolution induced by the 
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time-dependent Hamiltonian linearly interpolating between H{{j — 1)/ J) and H{j/ J) over a 
period of r/ J. Let Mj consist of backward, forward, and backward evolutions, namely. 



Mj = exp 



iH 



Uj exp 



iH 



J) 2J 



(7) 



J J 2J 

Our full state-preparation process is n/=i ^j- The dynamical phases converge to zero as J — )■ 
oo, while the adiabatic change of eigenbasis is undisturbed (§A.5). 

After the system has evolved for a period in which scattering occurs, measurement is per- 
formed as follows. The interaction is adiabatically turned off, through the time-reversed version 
of the turn-on described above. Once we return to the free theory, we can measure the number 
operators of the momentum modes, using the method of phase estimation, that is, by simulating 
giL <^aj,apt j-Qj. various values of t and Fourier transforming the results {31). 

Having described how, once discretized, a quantum field theory becomes essentially an 
ordinary many-body quantum-mechanical system, whose evolution can be efficiently simulated 
on quantum computers by combining established primitives, we now consider discretization 
errors. To analyse the errors introduced to our simulation by discretization, we use methods 
of effective field theory, a well-developed formalism underlying our modem understanding of 
quantum field theory. 

In its regime of validity, typically below a particular energy scale, an effective field theory 
(EFT) reproduces the behavior of the full (that is, fundamental) theory under consideration: it 
can be regarded as the low-energy limit of that theory. An EFT for a full theory is thus some- 
what analogous to a Taylor series for a function. It involves an expansion in some suitable small 
parameter, so that, although it consists of infinitely many terms, higher-order terms are increas- 
ingly suppressed. Thus, the series can be truncated, with corresponding finite and controllable 
errors. 

We apply this framework to analyse the effect of discretizing the spatial dimensions of the 
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continuum 0^ quantum field theory. The discretized Lagrangian can be thought of as the leading 
contribution (denoted by to an effective field theory. From the leading operators left out 
we can thus infer the scaling of the error associated with a non-zero lattice spacing, a. 

The full (untruncated) effective Lagrangian will have every coupling respecting the — > —(f) 
symmetry, and so will take the form 

£eff = C^'^ + + <^ct>'d'ct> + + • • • ■ (8) 

This can be simplified. First, the chain rule and integration by parts (with boundary terms 
dropped) can be used to write any operator with two derivatives acting on different fields in 
the form 0"c?^0. For example, 0^c?^09'*0 = |c?^(0^)c?''0 —^(p^d'^cf) . Such an operator 
can then be simplified via the equation of motion (32, 33). If this were the equation of motion 
of the continuum theory, any derivative operator would then be completely eliminated. In the 
discretized theory, however, the equation of motion is modified and there are residual, Lorentz- 
violating operators. In fact, because the difference operators in the discretized theory are only 
approximately equal to the derivatives in the continuum theory, the simplest Lorentz-violating 
operators are induced purely by discretization. 

In units where h — c — 1, all quantities have units of some power of mass. The mass 
dimensions (denoted by [.]) of the field and coupling in D — d+ 1 spacetime dimensions are 
[(f)] = ^ and [A] = 4 - D, which imply that 

[c] = 6 - 2D , [c"] = 8 -3D. (9) 

In D — A dimensions, [c] = —2 and [c"] — —4. Since the only relevant dimensionful parameter 
is the lattice spacing, that is, A ~ vr/a, this means that c ~ and c" ~ a*^. We see then that, 
of the operators not included in the Lagrangian (f>^ is more significant than 0^", for n > 3. 

In D — 2,3, the scaling of the coefficients with a is somewhat less obvious, because now 
the coupling A provides another dimensionful parameter. To obtain the scaling of c, one should 
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consider the Feynman diagram that generates the corresponding operator. This involves three 

(j)^ vertices, so 



(Other diagrams involve higher powers of A and hence are suppressed by higher powers of a.) 
Likewise, the coefficient of 0® will scale as X'^a^~^, which means that it is suppressed by 
relative to the coefficient of 0^. 

The effective field theory thus consists of three different classes of operators: operators of 
the form 0^", Lorentz-violating operators arising solely from discretization effects, and Lorentz- 
violating operators due to discretization and quantum effects. These are shown with the scaling 
of their coefficients in Table 2. At strong coupling, the operators and their scaling remain 
the same at the scale of the matching of the full theory on to the EFT, although the explicit 
coefficients are no longer calculable. However, the running of the coefficients down to lower 
energies is determined by their so-called anomalous dimensions, which depend on the coupling 
strength. These anomalous dimensions modify the scaling; at weak coupling the modification 
is small, but at strong coupling it could be larger. (Still, the scaling will remain polynomial.) 

From Table 2, one sees that the dominant discretization errors scale as a'^ in D — 2,3, 4. 
(In D = 2, 3, errors of type II dominate. In D = 4, errors of types I and II each scale as a^.) 
These error terms shift scattering probabilities by ±e where e = O(a^). Limiting e determines a 
and V — ^, which then determine the overall complexity of the quantum simulation algorithm 
described in Eq. 1 and Table 1. 

In conclusion, we have shown that quantum computers can efficiently calculate scattering 
probabilities in 0^ theory to arbitrary precision at both weak and strong coupling. Known 
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Class 


Operators 


Scaling of coupling 


1 


(n > 3) 




II 


0^2^ (/ > 2) 


a2'-2 


m 


0' > 1> ^ > 2) 


Xj+l(j2j+2l+2-D 



Table 2: Effective field theory operators fall into three classes. The general operator in each 
class is shown, with the canonical scaling of its coefficient in D spacetime dimensions. Here, 

classical algorithms take exponential time to do this in the strong-coupling and high-precision 
regimes. In addition to establishing a new exponential quantum speedup, our algorithm intro- 
duces several new techniques. These lead the way towards a quantum algorithm for simulating 
the Standard Model of particle physics, which has new features, such as chiral fermions and 
gauge interactions. Such an algorithm would establish that, except for quantum-gravity effects, 
the standard quantum circuit model suffices to capture completely the computational power of 
our universe. 
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A Supplementary Material 

A.l Steps of Algorithm and Comments 

Our quantum algorithm works by the following sequence of steps. 

1. Prepare the free vacuum. Improving upon the efficiency of earlier, more general, state- 
construction methods (6, 34), Kitaev and Webb developed a quantum algorithm for con- 
structing multivariate Gaussian superpositions (27). For large V, the dominant cost in 
Kitaev and Webb's method for producing V-dimensional multivariate Gaussians is the 
computation of the LDL^ decomposition of the inverse covariance matrix, where L is a 
unit lower triangular matrix, and D is a diagonal matrix. This can be done in 0(V^'^^^) 
time with established classical methods (28, 29). (The notation f{n) = 0{g{n)) means 
/(n) = 0{g{n) log'^(n)) for some constant c.) The computation of the matrix elements 
of the covariance matrix itself is easy because, for large V, the sum 



defining the propagator of the lattice theory is well approximated by an easily evaluated 
integral. 

2. Excite wavepackets. The span of |vac(0))|0) and |V')|1) is an invariant subspace, on 
which acts as 




(11) 



//v-|vac(0))|0) 



\m), 



(12) 



vac(0))|0) . 



(13) 



Thus, 



g-iH^-/2|^ac(0))|0) = -i|V')|l). 



(14) 
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Hence, by simulating a time evolution according to the Hamiltonian H^, we obtain the 
desired wavepacket state |^), up to an irrelevant global phase and extra qubit, which can 
be discarded. After rewriting in terms of the operators 0(x) and 7r(x), one sees that 
simulating is a very similar task to simulating H, and can be done with the same 
techniques. 

The only errors introduced at this step are due to the finite separation distance S be- 
tween wavepackets, and are of order e ~ e~^^'^. (However, our wavepackets have a 
constant spread in momentum, and thus differ from the idealization of particles with 
precisely defined momenta.) The wavepacket preparation thus has complexity scaling 
linearly with njn, the number of particles being prepared, and necessitates a dependence 
y-niJog(l/e). 

3. Adiabatically turn on the interaction. For < s < 1, let 

with Ao(0) = 0. Uj is the unitary time evolution induced by Hit/r) from t — ^-j io 
t = (2^, namely. 







Uj = rjexp 


—i / H(s)rds 




Jj/J 



where r{-} indicates the time-ordered product. We suppress the dynamical phases by 
choosing J to be sufficiently large. The choice of a suitable "path" Ao(s),mQ(s), and 
the complexity of this state-preparation process depends in a complicated manner on the 
parameters in H (§ A.5). 

4. Simulate Hamiltonian time evolution. 

5. Adiabatically turn off the interaction. The adiabatic turn-off of the coupling is simply 
the time-reversed version of the adiabatic turn-on. 
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7r(x)^ + l(V„0)^(x) + im2(.)0(x) 



(15) 



6. Measure occupation numbers of momentum modes. For a given p, measurement of 

L~'^aJ,ap by phase estimation can be implemented with O {v^^^^ quantum gates via 
a A;*^-order Suzuki-Trotter formula. Furthermore, if we instead simulate localized de- 
tectors, the computational cost becomes independent of V (much as the computational 
cost of creating local wavepackets is independent of V), but the momentum resolution 
becomes lower, as dictated by the uncertainty principle. 

The allowable rate of adiabatic increase of the coupling constant during state preparation 
is determined by the physical mass of the theory. In the weakly coupled case, this can be 
calculated perturbatively. In the strongly coupled case, such a calculation is no longer possible. 
Thus one is left with the problem of determining how fast one can perform the adiabatic state 
preparation without introducing errors. Fortunately, one can easily calculate the mass on a 
quantum computer, as follows. First, one adiabatically prepares the interacting vacuum state 
at some small Aq, and measures the energy of the vacuum using phase estimation. The speed 
at which to increase Aq can be chosen perturbatively for this small value of Aq. Next, one 
adiabatically prepares the state with a single zero-momentum particle at the same value of Aq, 
and measures its energy using phase estimation. Subtracting these values yields the physical 
mass. This value of the physical mass provides guidance as to the speed of adiabatic increase 
of the coupling to reach a slightly higher Aq. Repeating this process for successively higher 
Aq allows one to reach strong coupling, while always having an estimate of mass by which to 
choose a safe speed for adiabatic state preparation. In addition, mapping out the physical mass 
as a function of bare parameters (hence, for example, mapping out the phase diagram) may be 
of independent interest. 
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A.2 Efficiency 

To quantify the precision of a simulation, we demand that the probability of a given scattering 
event in the simulation differ from the true physical probability by no more than ±e. There are 
various sources of error: discretization of space, Trotter approximations, imperfect adiabaticity, 
discretization and cutoff of the field at each site, and imperfect spatial separation of particles 
in the asymptotic in and out states. In a theory with a non-zero mass, errors due to imperfect 
particle separation shrink exponentially with distance. Thus, V needs to scale only logarith- 
mically with e. Similarly, by the analysis of § A.4, the number of qubits per site scales only 
logarithmically with e. By Eq. 110, the errors resulting from use of a A;*^ -order Suzuki-Trotter 
formula with n timesteps are e ~ n"^''. Thus, the complexity scales as e"^/^'''. For large k, the 
dominant contributions to scaling with e are spatial discretization and imperfect adiabaticity. 

The effect of spatial discretization is captured by (infinitely many) additional terms in the 
effective Hamiltonian. Truncation of these terms alters the calculated probability of scattering 
events. In particular, the two dominant extra terms in the effective Hamiltonian are ^ ■ (f)df(f) 
and (f)^ terms, arising from discretization of (Va^)^ and quantum effects, respectively. The 
coefficient of the '^^(pdfcl) term is O(a^), and the coefficient of the 0^ term is 0{a^~'^), so 
that the former dominates for o? = 1,2, whereas the latter makes a comparable contribution for 
d — 3. Thus, the overall discretization error is 

e = 0(a^), d=l,2,3. (17) 

(To improve the scaling, one can use better finite differences to approximate the derivative, 
and/or include the (j)^ operator. However, renormalization and mixing of the coefficients make 
this idea more complicated than it is in standard numerical analysis.) 

The diabatic errors at weak coupling are estimated and summarized in § A.5.1. The errors 
are quantified by a probability e of observing stray particles. Substituting the a ~ depen- 
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l+<i/2+o(l) 

Gadiabatic ~ I " I , 1,2,3 (18) 



dence from Eq. 17 into Eq. 95 yields 

scaling for the adiabatic state preparation. We use little-o notation to convey precisely that we 
are neglecting both logarithmic factors and contributions to the exponent that become arbitrarily 
small as we use higher-order Suzuki-Trotter formulae. The other slow part of the algorithm is 
the preparation of the free vacuum. This scales as 

Gprep = 0(V2-3^^) = 0(a-2-376'^) = ©(e-i-iss^^) , (19) 

where the last equality follows from Eq. 17. Thus, ind—1 the adiabatic state preparation is the 

dominant cost, whereas in d = 2, 3 the preparation of the free vacuum dominates. This leaves a 
final asymptotic scaling of 

_ I 

(jtotal — 0(Gadiabatic + G'prep) — \ (-) 

i &' 

The number of quantum gates used to simulate the strongly coupled theory has scaling in 
^/{K — ^0) and p that is dominated by adiabatic state preparation (§ A.5.2). We also estimate 
scaling with riout as follows. For two incoming particles with momenta p and — p, the maxi- 
mum number of kinematically allowed outgoing particles is riout ~ P- For continuum behavior, 
p = 7] /a for constant 77 -C 1. Furthermore, one needs V ~ nout to obtain good asymptotic out 
states separated by a distance of at least ~ 1/mo. Thus, V ~ n^ut' so one needs 'i\ut^^'^^^^ 
gates to prepare the free vacuum and, by Eq. 107, n^^^^'^"^^'' gates to reach the interacting the- 
ory adiabatically. (The adiabatic turn-off takes no longer than the adiabatic turn-on.) Hence the 
total scaling in riout is dominated by preparation of the free vacuum in three-dimensional space- 
time, but by adiabatic turn-on in two-dimensional spacetime. These results are summarized in 
Table 1. 



1.5+0(1) 


d = 


1, 




2.376+o(l) 


d = 


2, 


(20) 


3.564+o(l) 

1 


d = 


3. 





'* Whether we use Eq. 95 or Eq. 96 affects only the scahng with V. 
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A.3 Mass Renormalization 

The physical, or renormalized, mass as a function of the coupling features prominently in our 
calculations. For the weak-coupling regime, its form is obtained by perturbation theory. For the 
strong-coupling regime, we use its known behavior near the phase transition. 

At first order in the coupling, the shift of the squared mass is given by i times the one-loop 
Feynman diagram 



The calculation of these diagrams is quite analogous to standard calculations in perturbative 
quantum field theory, but there are a couple of differences. First, the propagator is different 
because of the discretization. Secondly, integrals over components 1, . . . , (but not component 
0) of loop momenta are cut off by vr/a, that is, the lattice spacing acts as an ultraviolet regulator. 
These differences alter the nature the integrals and hence what methods can be used to evaluate 
them. 

The existence of a phase transition in the cp'*' theory in D = 2 or 3 spacetime dimensions 
was shown rigorously in (20-22). As the system approaches it, thermodynamic functions and 
correlation functions exhibit power-law behavior, as is characteristic of a second-order phase 
transition. In particular, for constant ml. 




(21) 



At second order, there is also a contribution from the two-loop diagram 




(22) 



m ~ 



u 



(23) 



where Ac, the critical value of the coupling, depends on ml 
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Empirically, it has been found that systems with second-order phase transitions can be clas- 
sified into universality classes. Within each class, critical exponents are universal, taking the 
same values for all systems. (This universality is explained by the concept of the renormaliza- 
tion group.) The 0^ theory is believed to be in the same universality class as the Ising model, 
for which 



The value above for D = 3 has also been obtained directly in the 0^ theory by Borel resumma- 
tion (35). 

In D = 4 dimensions, in contrast, the believed triviality of the continuum 0^ theory im- 
plies that there is no non-trivial fixed point of the renormalization group and hence no phase 
transition as one varies (ml, Aq). Moreover, triviality places bounds on the maximum value of 
the renormalized coupling (36). In particular, strong coupling requires pa to be in the 

continuum-like regime, renormalized perturbation theory should be valid. 

A.4 Representation by Qubits 

The required number of qubits per site is 



In this section we show that one can simulate processes at energy scale E, while maintain- 
ing 1 — e fidelity to the exact state, with jib logarithmic in 1/a, 1/e, and V. Our analysis is 
nonperturbative, and thus applies equally to strongly and weakly coupled 0^ theory. 
Let \tjj)he the state, expressed in the field representation, namely, 




L> = 2, 
D = 3. 



(24) 



^6 = log (1 214>^sjk/ Sct>\) . 



(25) 




(26) 



and let 




(27) 
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Then 

(^/^l^/'cut) = / #1 . . . / #v p{4>i, ■■■,<Pv), (28) 



rmax 



where p is the probabiUty distribution 

p(0i,...>v) = |V'(0i,--->v)r- (29) 

In other words, (V'lV'cut) = 1 — Pout> where pout is the probability that at least one of 0i, . . . , 0v 
is out of the range [— 0max, 4>max\- By the union bound (Pr(A U B) < Pic{A) + Pr(S)), 

(V'lV'cut) > 1 - V max Pout (x) , (30) 

where Pout(x) is the probability that 0(x) is out of the range [— 0max, 0max]- 

Let //<^(x) and (70(x) denote the mean and standard deviation of 0(x) determined by p. By 
Chebyshev's inequality, choosing ^^ax = A*<)!)(x) + c(7<^(x) ensures 

Pout(x) < \ . (31) 



C2 



Thus, choosing 



0max = O ( max I + y^a^(x) ) ) (32) 



ensures (V'lV'cut) > 1 — e- 

Next, we observe the following, which is straightforward to prove. 

Proposition 1 Let p and q be Hermitian operators on (M) obeying the canonical commuta- 
tion relation [p, q\ — i\. Then the eigenbasis ofp is the Fourier transform of the eigenbasis of 

By Proposition 1, the eigenbasis of a'^7r{x) is the Fourier transform of the eigenbasis of 
0(x). Thus, discretizing 0(x) in increments of ^^(x) is roughly equivalent to the truncation 

-TTmax < 7r(x) < TT^ax, whcrC 

TTmax T7 ■ (33) 

a'*^0(x) 
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By the same argument used to choose 0max, choosing 

TTmax = O ^max ^/X^(x) + (Tt,{^) ) ) (34) 

ensures fideUty 1 — e between and its truncated and discretized version. 

To obtain useful bounds on 0max and TTmax^ we must bound /x</,(x), C0(x)> A*7r(x)> and cr7r(x)- To 
this end, we make the following straightforward observation. 

Proposition2 Let M be a Hermitian operator and let be a quantum state. Then\{%l)\M\%l))\ < 

Proof: For brevity, let {Q) = (V'lQlV') for any observable Q. The operator (M - {M)lf is 
positive semidefinite. Thus, 

< {{M-{M)tf) (35) 
= (M^ -2(M)M+ (M)2l) (36) 
= (M^)-(M)^ (37) 



□ 



Applied to the definitions 



i^m = (^I^^WK/^): (38) 
<T^(x) = v^(v#WW^^WWK^, (39) 

/^Mx) = (^k(x)lii'). (40) 



(j,(x) = V(V'k(x)2|V^)-(V.|7r(x)|V')2, (41) 



Proposition 2 implies that //</,(x) and (T0(x) are each at most •>/ ('0|0(x)2|'0), and /i7r(x) and (T7r(x) 
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are each at most \J (■0|7r(x)2|'0). Thus, by Eq. 32 and Eq. 34, 



O I max -i / — (■?/'|7r(x)2| 



(42) 

(43) 



so that, by Eq. 25 and Eq. 33, 



rib = O log a'^- max V(^|7r(x)2|^/;)(^|0(y)2|^) 



(44) 



To establish logarithmic scaling of nt, we need only prove polynomial upper bounds on 
(V'|0(x)2|^) and (V'|7r(x)2|'0). Rather than making a physical estimate of these expectation 
values, we prove simple upper bounds that are probably quite loose. In the adiabatic state prepa- 
ration described in § A.5, the parameters mg and Aq are varied. The following two propositions 
cover all the combinations of parameters used in the adiabatic preparation and subsequent scat- 
tering of both strongly and weakly coupled wavepackets. 

Proposition 3 Let H be of the form shown in Eq. 15. Suppose itIq > and Xq > 0. Let be 

any state of the field such that {iIj\H\iIj) < E. Then Vx G VL, 

2E 



d^2 ' 







a" 



(45) 
(46) 



Proof: 



E > 



\H\ 



-7r(x)2 + ^(V„0)2(x) + ^0(x)2 + ^0(x 



An 



4! 



> (V^|a'^^0(x)», 



(47) 
(48) 

(49) 
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where the last inequaUty follows because all of the operators we have dropped are positive 
semidefinite. This establishes Eq. 45. Similarly, we can drop all but the 7r(x) term from the 
right-hand side of Eq. 48, leaving 

£^>(V'|a''^7r(x)2|V^), (50) 

which establishes Eq. 46. □ 

Proposition 4 Let H be of the form shown in Eq. 15. Suppose ml < and Aq > 0. Let be 
any state of the field such that < E. Then Vx e Q, 



(^|0(x)2|^) < [e+ ^ ° h (51) 



Ao y A^ Xoa<^ V 2Ao 

nd 



{i^H^m < + , (52) 



where V is the physical volume. 
Proof: The operator 

t/(x) = ^<Pi^r + ^0(x)^ (53) 

is sufficiently simple that we can directly calculate its minimal eigenvalue C/min- lfml<0 and 
A > 0, then 

C/„.in = ■ (54) 



Thus, for any state l'^) 

2Ao 



{^\J2a'Ui^M>^^. (55) 
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Hence, recalling Eq. 15, we obtain 



E > i^m) 



An 



4! 



> 



-TT X 



^7r(xf + i(V„0)2(x) 



2Ao 



(56) 
(57) 

(58) 

(59) 



2 2Ao 

Eq. 58 follows from Eq. 55. Eq. 59 holds (for any choice of x) because all of the operators we 
have dropped are positive semidefinite. This establishes Eq. 52. 

Similarly, dropping positive operators from Eq. 57 and using Eq. 55 yield, for any x. 



Applying Proposition 2 with M = 0(x)^ shows that {ijj\(j){x.)'^\ijj) > {ijj\4>{x)'^\il!y. Thus, 



(60) 



Via the quadratic formula, this implies Eq. 51. 



< E + 



3{V - a'^)m^ 
2A^ 



(61) 

□ 



A.5 Adiabatic Preparation of Interacting Wavepackets 



In this section, we analyze the adiabatic state-preparation procedure. To analyze the error due to 
finite T and J, we consider the process of preparing a single-particle wavepacket. The procedure 
performs similarly in preparing wavepackets for multiple particles provided the particles are 
separated by more than the characteristic length 1/m of the interaction. 

The phase induced by Mj on the momentum-p eigenstate of H{s) (with energy Ep{s)) is 



J + 1 
J 



dsEp{s) . 



(62) 
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Taylor expanding Ep about s — {j + ^)/J yields 



+ o{j-'). 



12J3 ds^ 



(63) 



Thus the total phase induced is 



j-i 

j=0 



12J2 



r dEr, 



(64) 

(65) 
(66) 



12J2 ds 

where the approximation holds for large J. For a Lorentz-invariant theory, Ep{s) must take the 
form 

Ep{s) = + (67) 

This should be a good approximation for the lattice theory provided the particle momentum 
satisfies p <^ 1/a. Substituting Eq. 67 into Eq. 66 yields 



0{P) 



ds 



(68) 



24J2 ^p2 + m^(^s) 

Next, we consider the effect of this phase shift on a wavepacket centered around momentum 
p. If the wavepacket is narrowly concentrated in momentum, then we can Taylor expand 9{p) 
to first order about p: 

e{p)c^0{p)+V-ip-p), (69) 



where 



V- — 

dp 



(70) 



The phase shift e*^ induces a translation (in position space) of any wavepacket by a dis- 
tance v. (The second-order term in the Taylor expansion induces broadening.) From Eq. 70 
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iTT-o symmetric , 
phase B / 




broken-symmetry 



^ A 



phase 



Figure 1 : The dashed line illustrates schematically the location of a quantum phase transition of 
0^ theory in two and three spacetime dimensions. A and B denote weakly and strongly coupled 
continuum-like theories, respectively. We prepare them adiabatically by following the arrows 
starting from the massive free theory (ttIq > 0, Aq = 0). To maintain adiabaticity, the path must 
not cross the quantum phase transition. 



and Eq. 68, we have 



We next determine the complexity by demanding that the propagation length V be restricted 
to some small constant, and that the probability of diabatic particle creation be small. Together, 
these criteria determine J and r. We can obtain a tighter bound in the perturbative case than in 
the general case, so we treat these separately. 

A.5.1 Weak Coupling 

In the perturbative continuum limit a — )■ 0, ml is negative. For fixed small a, we can adiabati- 
cally approach a perturbative continuum-like theory by taking the straight-line path depicted in 
Fig. 1, namely, the following parameterization of Eq. 15: 



r\p\ ¥ 



V 



24J2 (p2+^2(5))3/2 



(71) 



s=0 




(m«)' + sAo/i 



Ao(s) 



sAo . 



(72) 
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Using perturbation theory (see diagram 21), one finds that it is particularly efficient to choose 



' logf^) +••• , (i = l, 
+ •••, d = 2, 

d = 3 , 



(2) 
^0 1 

167r^ a 

.(3) 



(73) 



_Zo_J_ J 

SO that, at first order in Aq, the physical mass remains fixed at mS^^ for all s. Here, Tq^^ 



25.379 ... and r, 



(3) 



112.948 In the perturbative regime, this should ensure that the path 



does not cross the quantum phase transition. 

To calculate the variation of physical mass with s, we must go to second order in Aq (see 
diagram 22). The result is 



where 



0(Ag/(mW)2) ^ ^^1^ 
m; = < O (Aglog(m(^)a)) , ci = 2 , 
0(A2/a2), d = 3. 



(74) 



(75) 



Substituting Eq. 74 into Eq. 7 1 yields 

rlpl 



<V. 



(76) 



12J'(p2 + (^(l))2 + ^2)3/2 

If we are considering a fixed physical process and using successively smaller a to achieve higher 
precision then, by Eq. 75, it suffices to choose J to scale as 



O 

6 
6 



AoI5 



gr \ 



(i = 2, 
d = 3. 



(77) 



Note that, for d = 3, J is suppressed by v^. This is because, as s increases, the (uncancelled) 
two-loop contribution to the physical mass makes the particle very heavy until s is very close to 
one. Hence, the particle propagates slowly, and less backward evolution is required. 
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To determine r, we next consider adiabaticity. Let H{s) be any Hamiltonian differentiable 
with respect to s. Let be an eigenstate H{s)\(f)i{s)) = Ei{s) separated by a non-zero 

energy gap for all s. Let {il^iit)) be the state obtained by Schrodinger time evolution according 
to H{t/T) with initial condition — The diabatic transition amplitude to any 

other eigenstate H{s)\(j)k{s)) = Ek{s)\(f)k{s)) (k ^ I) is (57) 

iMsMirs)) ~ rrf^ ^y;^)'^ '^'(^)) e-(^^W-^'(-)) (1 + 0(l/r)) . (78) 
Jo Ei{a) - hk[a) 

(The integrand is made well-defined by the phase convention = 0.) Here, 

ipi{s)= f daEi{a). (79) 
In the case that Ei, E^, and (0fe|^|0z) are s-independent, this integral gives 

(MsMirs)) ~ (1 - e^^^^-^')^) ^^^^^(1 + Oil/r')) . (80) 

In the case that these quantities are approximately s-independent, Eq. 80 should hold as an 
approximation. 

In reality, we wish to prepare a wavepacket state, not an eigenstate. However, the wavepacket 
is well separated from other particles and narrowly concentrated in momentum space. Thus, we 
shall approximate it as an eigenstate |0;(s)). Furthermore, by our choice of path, the energy gap 
is kept constant to first order in the coupling, and thus Eq. 80 should be a good approximation 
to Eq. 78. 

Summing the transition amplitudes to some state from the J steps in in our preparation 
process, and applying the triangle inequality^ yield the following: 



{Mj/j)mMj/j)) 



{E,{j/J) - E,{3/J)Y 



(81) 



^The 0{J) scaling obtained by the triangle inequality can be confirmed by a more detailed calculation taking 
into account the relative phases of the contributions to the total transition amplitude. 
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The j — term in this sum can be evaluated exactly, because it arises from the free the- 
ory. At j 7^ the theory is no longer exactly solvable. However, one obtains the lowest-order 
contribution to the matrix element (pi, p2, Ps, P4; s — l|0^|vac(l)) in renormalized perturba- 
tion theory simply by taking the j — expression and replacing mo with the physical mass 
and Ao with the physical coupling. Our adiabatic path Eq. 72 is designed so that the physical 
mass at s = 1 matches the bare mass at j (at least to first order in Aq). Furthermore, the 
physical coupling differs from the bare coupling only by a logarithmically divergent (in a) cor- 
rection for (i = 3 and non-divergent corrections for d = 1, 2. Thus we can make the following 
approximation: 



Diabatic errors come in two types, creation of particles from the vacuum, and splitting of 
the incoming particles. The matrix element in the numerator of Eq. 82 can correspondingly 
be decomposed as the sum of two contributions. We first consider particle creation from the 
vacuum, approximating \ (j)j{s)) as |vac(s)). 



Substituting this into the numerator of Eq. 82, setting \(f)i{0)) — |vac(0)), and expanding in 
terms of creation and annihilation operators show that the only potentially non-zero transition 
amplitudes are to states |0fc(O)) of two or four particles. The transition amplitude to states of 
four particles arise solely from the 0^ term in The transition amplitude to states of two 
particles has contributions from the 0^ term and the 0^ term in These actually cancel, 
because of our choice of /x. (Note that this requires tuning of //.) At s = 0, the numerator of 
Eq. 82 is therefore the following: 




(82) 



By Eq. 72, 




(83) 



(Pi,P2,P3,P4|7^XlaV(x)|vac(0)) 




(84) 
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We obtain the probability of excitation due to creation of four particles from the vacuum by 

squaring the amplitude estimated above, and then summing over all allowed combinations of 
the four outgoing momenta: 

^0"P1+P2+P3+P4,0 



create 



y 

PI pf^P4&r ^^'^^('^(Pi) + '^(P2) + ^{Ps) + u;(p4))^a;(pi)a;(p2)a;(p3)a;(p4) 



(85) 

This sum is difficult to evaluate exactly; instead, we shall simply estimate its asymptotic scaling. 
The question is, with which parameter should we consider scaling? There are at least three 
regimes in which classical methods for computing scattering amplitudes break down or are 
inefficient: strong coupling, large numbers of external particles, and high precision. In this 
section we are considering only weak coupling (that is, A/m^~^ ^ 1), leaving discussion of 
strong coupling until the next section. For an asymptotically large number of external particles, 
the efficiency of our algorithm depends upon strong coupling, for the following reason. A 
connected Feynman diagram involving n external particles must have at least v — 0{n) vertices, 
so the amplitude for such a process is suppressed by a factor of ( ^4-0 )^, where E is the energy 
scale of the process. Since E > m, many-particle scattering events are exponentially rare at 
weak coupling, and thus cannot be efficiently observed in experiments or simulations. This 
leaves the high-precision frontier. Recall that the perturbation series used in quantum field 
theory are asymptotic but not convergent. Thus, perturbative methods cannot be extended to 
arbitrarily high precision. 

Hence, in this section we consider the quantum gate complexity of achieving arbitrarily high 
precision. To do so, one chooses a small to obtain small discretization errors, V large to obtain 
better particle separation, r long to improve adiabaticity, and J large enough to limit unwanted 
particle propagation as the interaction is turned on. Thus, we wish to know the scaling of Pcreate 
with a, T, V, and J. In this context, we consider m. A, and |pi| to be constants. 
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We now estimate the scaling of Pcreate as a — > 0. 

J2 

^create ~ y2^2 5Z 

pi,P2,P3er 



^0 



a;(pi) + uj{p2) + i^iPs) + '^(-Pi - P2 - P3))^^(Pi)^(P2)'^(P3)'^(-Pi - P2 - Pa) 
3J2 



pi>P2,p3Gr 

|Pl|>|P2|,|P3| 



\2 



(a;(pi) + w(p2) + a;(p3) + a;(-pi - P2 - P3))''w(pi)a;(p2)a;(p3)c^(-pi - P2 - P3) 

^ n .r ^(Pl)'^(P2)^(P3) 

|Pl|>|P2|,|P3| 
72 \2 



Jr ir ir tj(pi)t>cu(p2)cj(p3) 
O(^) , d=l,2, 

By Eq. 77 and Eq. 86, 

Pcreate= O (^) , d = 1, 2, 3 . (87) 

Next, we consider the process in which the time dependence of the 0^ term causes a single 
particle to split into three. For this process, the relevant matrix element is 

/ I "^0 dlA/ M \ -^0<^P2+P3+P4,Pl /00\ 

P2,P3,P4 hrZ^Q X Pi ^ ^ ^ ^ ==, (88) 

4! 4\/Vc^(pi)c^(p2)c^(p3)c^(p4) 

where pi is the momentum of the incoming particle. By our choice of path, the physical mass 
is s-independent to first order in the coupling, and the s dependence of the coupling is only 
logarithmically divergent as a — > 0. Thus, by Eq. 81, 

p j.^ ^0"P2+P3+P4,Pl ^gQ-j 

spirt ^2y2 (a;(p2) + cu(p3) + Cc;(p4) - u[Y>\)f^[V\)^[V2)^{v-i)^{'?^ ' 

P2,P3,P4fci 
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(86) 



Let us now examine the divergence structure of Pgpiit as a —> 0. In the Umit of large volume, 

the sum converges to the following integral: 

— [ d'^P2 [ d'^P3 — 

t2 7r Jt (^(P2) + ^(Pa) + w(Pi - P2 - Pa) - w(Pi))^^^(Pi)w(P2)w(P3)w(Pi - P2 - P3) 

(90) 

If this were divergent as a — > 0, then by approximating the integrand with its value at large |p2| 
and I Pal, we would be able to isolate the divergence: 

Pspiit ^ -4^, f d'p2 [ d% J . (91) 

T^w(pi) 7r Jr (IP2I + IPsI + |P2 +P3|)^|P2||P3||P2 + P3I 

However, for d = 1, 2, 3 this is convergent as a — >^ 0. Thus, recalling Eq. 77, we obtain 

^spU-U^^^J -| o{^) , d^3. ^^^^ 

We can consider two criteria regarding diabatic particle creation. If our detectors are local- 
ized, we may be able to tolerate a low constant density of stray particles created during state 
preparation. This background is similar to that encountered in experiments, and may not in- 
validate conclusions from the simulation. Alternatively, one could adopt a strict criterion by 
demanding that, with high probability, not even one stray particle is created in the volume being 
simulated during state preparation. This strict criterion can be quantified by demanding that the 
adiabatically produced state has an inner product of at least 1 — e with the exact state. This 
parameter e is thus directly comparable with that used in § A.4, and the two sources of error can 
be added. Applying the strict criterion, we demand that Pgpiit and Pcreate each be of order e, and 
obtain 

Tsw = O (^^^ , d= 1,2,3. (93) 
Applying the more lenient criterion that Pcreate/^ and Pgpiit each be of order e yields 

Tlenient ^ O ( -) , d ^ 1,2,3 . (94) 
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For a /c**^ -order Suzuki-Trotter formula, the asymptotic scaling of the total number of gates 
needed for adiabatic state preparation is O ( (Vr)^+^ ) = O ( (V'r/a'^)^+^) ) . Thus, 



A.5.2 Strong Coupling 

In two and three spacetime dimensions, we can obtain a strongly coupled (that is, nonperturba- 
tive) field theory by approaching the phase transition (§ A. 3). As in the case of weak coupling, 
the necessary time for adiabatic state preparation depends on various physical parameters of 
the system being simulated, including the momentum of the incoming particles, the volume, 
the strength of the final coupling, the number of spatial dimensions, and the physical mass. To 
keep the discussion concise, we restrict our discussion to the case of ultrarelativistic incoming 
particles, with coupling strength close to the critical value. Under these conditions, the incom- 
ing particles can produce a shower of many (uout ~ p/m) outgoing particles. Because of the 
strong coupling, perturbation theory is inapplicable, and, even if it could be used, would take 
exponential computation in the number of outgoing particles. 

In the strongly coupled case, we vary the Hamiltonian 15 with s by keeping the bare mass 
constant at mo and setting the bare coupling to sXq. We choose Aq only slightly below the critical 
value Ac, so that at s = 1 the system closely approaches the phase transition, as illustrated in 
Fig. 1. Examining Eq. 68 suggests that we can estimate phase errors by understanding the 
behavior of m^(s) at s = and s — 1, without needing to know exactly what happens in 
between. From Eq. 73, 




(95) 



(96) 




(97) 
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and, from Eq. 23 and Eq. 24, 



ds 

Thus, Eq. 71 yields 



s=l 



-2(A,-Ao) d=l, 
-1.26(Ae- Ao)°-26 d = 2. 



0.26 J _ o (98) 



under the assumption that (Ac — Aq) is very small. 

The result 68 rests on two approximations, a Taylor expansion to second order in Eq. 63, 
and an approximation of a sum by an integral in Eq. 66. The validity conditions for these 
approximations become most stringent at s = 1, where the derivatives of with respect to s 
become large. Working out the 0( J term in Eq. 68 at s = 1, one finds that it will be much 
smaller than the 0( J~^) term at s = 1 provided 

J » . (100) 

Similarly, higher-order terms in the Taylor expansion are suppressed by additional powers of 
j(A -Ao) • criterion 100 also suffices to justify the approximation of the sum by an integral 
in Eq. 66. 

We must next consider adiabaticity to determine r. In the ultrarelativistic limit, the relevant 



2 

energy gap 7 is ~ This takes its minimum value at s = 1, namely. 



d = 1 
d = 2. 



p 

Unlike in the perturbative case, we cannot make a detailed quantitative analysis, but under the 
condition 100, we should again be able to apply the traditional adiabatic criterion and obtain a 
diabatic transition amplitude scaling as Thus, to keep the error probability at some small 



constant e, we have 



T ~ ^ . (102) 
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We now consider asymptotic scaling with p for fixed Aq. To achieve continuum-Uke behavior 

we need a <^ ^- Thus Eq. 99 yields 

J ~ r^/'^pid-^)/'^ ^ d = 1, 2 . (103) 
Substituting Eq. 100 and Eq. 101 into Eq. 102, we see that we need 

T>p\ d=l,2. (104) 
Substituting Eq. 103 and Eq. 101 into Eq. 102, we see that we need 

r>p"+\ d=l,2. (105) 

The scaling r = 0{p'^^^) for d = 1,2 suffices to satisfy both conditions 104 and 105. Thus, by 
§A.6, the total number of gates scales as 

Strong - 0((yr)l+°(^y+^+«(l)) (106) 
= O (V^+oil)p2d+2+oil)^ ^ 

for d = 1,2. 

Next, we consider asymptotic scaling with (Ac — Aq) for fixed p. The J scaling as ^/T 
in Eq. 99 automatically satisfies the condition 100. Thus, we substitute Eq. 99 into Eq. 102, 
obtaining 

r~{ (108) 
Thus, using a k^^-order Suzuki-Trotter formula, we obtain 



Note that one could improve this scaling by choosing a more optimized adiabatic state-preparation 
schedule, which slows down as the gap gets smaller. 
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A.6 Suzuki-Trotter Formulae for Large Lattices 

It appears that, while scaUng with t has been thoroughly studied, little attention has been given 
to scaling of quantum simulation algorithms with the number of lattice sites V. Using a result 
of Suzuki and elementary Lie algebra theory, we derive linear scaling provided the Hamiltonian 
is local. 

For any even k and any pair of Hamiltonians A,B, 

^^iAait/n^iB/3it/n^iAa2t/n^il32Bt/n ^iAart/n>^^ _ ^i{A+B)t _|_ Q^^2k+1 y^2fe^ (110) 

where r — 1 + 5*^/^"^ and ai, . . . , a^, Pi, ■ ■ ■ , Pr-i are specially chosen coefficients such 
that ^^=1 otj = 1 and X]j=i l^j = ^ C^'^)- Thus, using the A;*'^-order Suzuki-Trotter formula 

/ 2fc+l \ 

(Eq. 110), one can simulate evolution for time t with O I i~2fc~ j quantum gates (4). To deter- 
mine the V scaling, we use the following standard theorem (cf. the Baker-Campbell-Hausdorff 
formula). 

Theorem 1 Let A and B be elements of a Lie algebra defined over any field of characteristic 
0. Then e^^e^ = , where C is a formal infinite sum of elements of the Lie algebra generated 
by A and B. 

A and B generate a Lie algebra by commutation and linear combination. Thus, without requir- 
ing any explicit calculation. Theorem 1 together with Eq. 110 implies 

(^^iA5,t/n^iB5,t/n ^iASrt/ny ^ ^i(A+B)t ^ /^^^^^f^k+l ^^2k ^ 0(^-(2fc+l)) ^ (J H) 

where A2jfc+i is a linear combination of nested commutators. In general, || A2fe+i|| could be as 
large as (max{||74||, However, by the canonical commutation relations, one sees 

that, for the pair of local Hamiltonians H^, if^, || A2fe+i|| = 0(V), for any fixed k. Thus, one 
needs only n — O (t^^V"^^ . Recalling the 0{V) cost for simulating each e^^'i'^^ or e^^'"^*', 
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one sees that the total number of gates scales as O ^(^V)^''"^^ . Note that this conclusion may 
be of general interest, as it applies to any lattice Hamiltonian for which non-neighboring terms 
commute. 

In the case of strong coupling, we care not only about how the number of gates scales with 
V, but also about scaling with p. In the presence of high-energy incoming particles, the field 
can have large distortions from its vacuum state. For example, if (V'|0(x) |^) is large, then local 
terms in /^2k+\W) such as 7r(x)0(x)^|z/)) can become large. We can obtain a heuristic upper 
bound on this effect by noting that, in the strongly coupled case, > 0, so each local term 
in if is a positive operator. Thus, if {'^\H\'^) < E, then the expectation value of each of the 
local terms is bounded above by E. Using £^ as a simple estimate of the maximum magnitude 
of a local term, we see that A2k+i |^), which is a sum of 0(V) terms, each of which is of degree 
2/i;-|- 1 in the local terms of H, has magnitude at most 0(V-E^*^"'"^), or in other words 0{Vp^''~^^). 
Recalling that a scales as a small multiple of 1/p, we see that A2fc+i|V') = 0{Vp'^'''^^'^'''). Thus, 
n = (9(^j)i+(i+'^)/2fc^i+i/2fc^ Each timestep requires 0(V) = 0{Vp'^) gates to implement. Thus, 
the overall scaling is 0{p'''~^^~^"'^^\tvy~^°^^^) quantum gates to simulate the strongly coupled 
theory at large p. 
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